import LinearAlgebra

function expandedCheb(ngp)
    res = zeros(ngp)
    if ngp == 1
        res[1] = 0.0
    else
        for k = 1:ngp
            res[k] = cos((2 * (ngp - k + 1) - 1.0) * pi / (2 * ngp)) / cos(pi / (2 * ngp))
        end
    end
    return res
end

function dqcoeff!(xlist::Vector{Float64}, maxorder::Int64, c::Array{Float64,3})
    #xlist: input, one-dimensional real list, with length n
    #maxorder : input, interger
    #c : output, differential quadrature coefficient array, with dimensions (maxorder+1, n, n)
    

    n = size(xlist, 1)
    c .= 0.0
    for i = 1:n
        for j = 1:n
            if i != j
                c[1, i, j] = Mp(xlist, i) / ((xlist[i] - xlist[j]) * Mp(xlist, j))
            end
        end
    end

    for i = 1:n
        c[1, i, i] = -sum(c[1, i, 1:i-1]) - sum(c[1, i, i+1:n])
    end

    for m = 2:maxorder
        for i = 1:n
            for j = 1:n
                if i != j
                    c[m, i, j] = m * (c[m-1, i, i] * c[1, i, j] - c[m-1, i, j] / (xlist[i] - xlist[j]))
                end
            end
        end
        for i = 1:n
            c[m, i, i] = -sum(c[m, i, 1:i-1]) - sum(c[m, i, i+1:n])
        end
    end

    for i = 1:n
        c[maxorder+1, i, i] = 1.0
    end
    return nothing
end

function Mp(xlist, i)
    #implicit none
    #real(kind=iwp) :: Mp, xlist(:)
    #integer(4) :: n, i, j
    n = size(xlist, 1)
    r = 1.0
    for j = 1:(i-1)
        r = r * (xlist[i] - xlist[j])
    end
    for j = (i+1):n
        r = r * (xlist[i] - xlist[j])
    end
    return r
end

function main()
    n = 30
    maxorder = 2
    x = expandedCheb(n)
    c = zeros(maxorder + 1, n, n)
    dqcoeff!(x, maxorder, c)

    y = sin.(x)

    yp = c[2,:,:] * y

    y1 = -sin.(x)

    diff = yp - y1
    print(LinearAlgebra.norm(diff))
end

main()